Alpha and Beta Diversity Analysis
Yixuan Yang
2024-12-29
Set Up
rm(list = ls())
library(tidyverse)
library(here)
library(phyloseq)
library(vegan)
library(rstatix)
set.seed(2024)
theme_set(theme_bw())
max.core <- parallel::detectCores()
ps.rare <- readRDS(here('data','following_study','ps_rarefied.rds'))
sample_data(ps.rare)$Shannon <- estimate_richness(ps.rare)$Shannon
# transform data into proportion
ps.rare.prop <- ps.rare %>% transform_sample_counts(function(x) x/sum(x))
sam <- data.frame(sample_data(ps.rare))Functions
plot_ord draws ordination plot for different factors
using plot_ordination function in phyloseq
package.
permanova performs permutational multivariate analysis
of variance (PERMANOVA) based on adonis2 function in
vegan package.
permdisp performs permutational analysis of multivariate
dispersions (PERMDISP) based on betadisper function in
vegan package.
plot_ord <- function(data, factor, method, distance){
data.ord <- ordinate(data, method = method, distance = distance)
p <- plot_ordination(data, data.ord, color = factor)
p <- p + stat_ellipse(type = "t",geom = "polygon",alpha = 0)
p <- p + ggtitle(str_c(factor,method,distance, sep = ' - '))
print(p)
}
permanova <- function(data, formula, method, permutations=1e4, strata = NULL, core = max.core){
message('PERMANOVA Model: ', method, '~', formula, '; Strata: ', ifelse(is_null(strata), 'None', as.character(strata)))
dist.matrix <- phyloseq::distance(data, method=method)
df <- data.frame(sample_data(data))
model <- as.formula(paste0('dist.matrix~', formula))
if (!is_null(strata)) {strata <- df[,strata]}
result <- adonis2(model,
data = df,
permutations=permutations,
strata = strata,
parallel = core,
by = 'term',
na.action = na.omit)
return(result)
}
permdisp <- function(data, group, method, permutations=1e4, pairwise = FALSE, core = max.core){
message('PERMDISP Model: ', method, '~', group)
dist.matrix <- phyloseq::distance(data, method=method)
df <- data.frame(sample_data(data))
beta.disp <- betadisper(dist.matrix, group = df[,group])
result <- permutest(beta.disp, permutations = permutations, pairwise = pairwise, type = 'centroid')
return(result)
}Start
In this section, we want to estimate the effect of different factors
on the microbial diversity. The factors we are focusing on are
Household, Epileptic.or.Control,
Breed.Group..1., Pheno.Y.N, Sex,
and Age..months.. We compare the species richness (Shannon
index) among different factors using ANOVA, and compare the centroid of
dissimilarity of microbial community between different groups using
PERMANOVA using the Bray-Curtis and weighted Unifrac distance, and
visualized using multi-dimensional scaling. PERMDISP was used to test
the homogeneity of multivariate dispersions among groups.
Household Effect
Alpha Diversity
ggplot(sam,aes(x = as.numeric(Household), y = Shannon, group = Household)) +
geom_point() + geom_line() + xlab('Household')
anova_test(Shannon~Household, data = sam, type = 1)## ANOVA Table (type I tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Household 48 49 1.591 0.054 0.609Here we see the Shannon diversity index is significantly different among households.
Beta Diversity
Bray-Curtis distance
permanova(ps.rare.prop, 'Household', 'bray')## PERMANOVA Model: bray~Household; Strata: None## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2401 0.68965 2.2685 9.999e-05 ***
## Residual 49 5.5082 0.31035
## Total 97 17.7482 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Weighted-Unifrac distance
permanova(ps.rare.prop, 'Household', 'wunifrac')## PERMANOVA Model: wunifrac~Household; Strata: None## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 1.90061 0.67646 2.1344 9.999e-05 ***
## Residual 49 0.90904 0.32354
## Total 97 2.80965 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Epileptic Effect
Alpha Diversity
ggplot(sam, aes(x = Epileptic.or.Control, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
anova_test(Shannon~Household+Epileptic.or.Control, data = sam, type = 1)## ANOVA Table (type I tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Household 48 48 1.560 0.063 0.609
## 2 Epileptic.or.Control 1 48 0.053 0.819 0.001Beta Diversity
Bray-Curtis distance
plot_ord(ps.rare.prop, 'Epileptic.or.Control','MDS','bray')
plot_ord(ps.rare.prop, 'Epileptic.or.Control','NMDS','bray')## Run 0 stress 0.2088385
## Run 1 stress 0.2038739
## ... New best solution
## ... Procrustes: rmse 0.05493282 max resid 0.4067546
## Run 2 stress 0.2215044
## Run 3 stress 0.2196132
## Run 4 stress 0.2048226
## Run 5 stress 0.2358109
## Run 6 stress 0.2302608
## Run 7 stress 0.2044114
## Run 8 stress 0.22333
## Run 9 stress 0.2040503
## ... Procrustes: rmse 0.03678772 max resid 0.3043432
## Run 10 stress 0.2209504
## Run 11 stress 0.2215212
## Run 12 stress 0.2115265
## Run 13 stress 0.2102944
## Run 14 stress 0.2093072
## Run 15 stress 0.2165779
## Run 16 stress 0.208752
## Run 17 stress 0.2215795
## Run 18 stress 0.2048226
## Run 19 stress 0.2040883
## ... Procrustes: rmse 0.0150325 max resid 0.09310887
## Run 20 stress 0.2200547
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 3: no. of iterations >= maxit
## 17: stress ratio > sratmax
permanova(ps.rare.prop, 'Epileptic.or.Control', 'bray', strata = 'Household')## PERMANOVA Model: bray~Epileptic.or.Control; Strata: Household## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Epileptic.or.Control 1 0.1585 0.00893 0.8653 0.09059 .
## Residual 96 17.5897 0.99107
## Total 97 17.7482 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1permdisp(ps.rare.prop, 'Epileptic.or.Control', 'bray')## PERMDISP Model: bray~Epileptic.or.Control##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00167 0.0016744 0.1118 10000 0.7439
## Residuals 96 1.43792 0.0149784Weighted-Unifrac distance
plot_ord(ps.rare.prop, 'Epileptic.or.Control','MDS','wunifrac')
plot_ord(ps.rare.prop, 'Epileptic.or.Control','NMDS','wunifrac')## Run 0 stress 0.1672646
## Run 1 stress 0.1719813
## Run 2 stress 0.1660392
## ... New best solution
## ... Procrustes: rmse 0.04624571 max resid 0.145459
## Run 3 stress 0.1676809
## Run 4 stress 0.1654881
## ... New best solution
## ... Procrustes: rmse 0.06177539 max resid 0.3530994
## Run 5 stress 0.1703517
## Run 6 stress 0.1695267
## Run 7 stress 0.1697312
## Run 8 stress 0.1635417
## ... New best solution
## ... Procrustes: rmse 0.0166237 max resid 0.1104734
## Run 9 stress 0.1767302
## Run 10 stress 0.1717797
## Run 11 stress 0.1739713
## Run 12 stress 0.1710472
## Run 13 stress 0.1733658
## Run 14 stress 0.1678347
## Run 15 stress 0.1691009
## Run 16 stress 0.1673507
## Run 17 stress 0.171604
## Run 18 stress 0.1730282
## Run 19 stress 0.1743485
## Run 20 stress 0.1679275
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 19: stress ratio > sratmax
permanova(ps.rare.prop, 'Epileptic.or.Control', 'wunifrac', strata = 'Household')## PERMANOVA Model: wunifrac~Epileptic.or.Control; Strata: Household## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Epileptic.or.Control 1 0.04072 0.01449 1.4118 0.0408 *
## Residual 96 2.76893 0.98551
## Total 97 2.80965 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1permdisp(ps.rare.prop, 'Epileptic.or.Control', 'wunifrac')## PERMDISP Model: wunifrac~Epileptic.or.Control##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00697 0.0069689 1.8477 10000 0.1768
## Residuals 96 0.36208 0.0037716Breed Effect
Alpha Diversity
sam.breed <- sam %>% filter(is.na(Breed.Group..1.) == FALSE)
ggplot(sam.breed) +
geom_point(aes(x = Breed.Group..1., y = Shannon, colour = Breed.Group..1.)) +
facet_wrap(~Epileptic.or.Control) +
theme(axis.text.x = element_blank(), axis.ticks.x.bottom = element_blank())
anova_test(Shannon~Household + Breed.Group..1., data = sam.breed, type = 1)## ANOVA Table (type I tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Household 45 39 1.576 0.075 0.645
## 2 Breed.Group..1. 4 39 2.930 0.033 * 0.231Beta Diversity
Bray-Curtis distance
plot_ord(ps.rare.prop, 'Breed.Group..1.','MDS','bray')## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.rare.prop, 'Breed.Group..1.','NMDS','bray')## Run 0 stress 0.2088385
## Run 1 stress 0.2247104
## Run 2 stress 0.2129973
## Run 3 stress 0.2097522
## Run 4 stress 0.2031357
## ... New best solution
## ... Procrustes: rmse 0.04265187 max resid 0.2822212
## Run 5 stress 0.2240573
## Run 6 stress 0.2056011
## Run 7 stress 0.2217547
## Run 8 stress 0.232074
## Run 9 stress 0.2039415
## Run 10 stress 0.2052108
## Run 11 stress 0.2029263
## ... New best solution
## ... Procrustes: rmse 0.0138755 max resid 0.08701649
## Run 12 stress 0.2192778
## Run 13 stress 0.2095446
## Run 14 stress 0.2060592
## Run 15 stress 0.2280706
## Run 16 stress 0.2312938
## Run 17 stress 0.2093494
## Run 18 stress 0.2221842
## Run 19 stress 0.2139507
## Run 20 stress 0.2029124
## ... New best solution
## ... Procrustes: rmse 0.002889526 max resid 0.02534544
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 6: no. of iterations >= maxit
## 13: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrmin## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
permanova(ps.rare.prop, 'Household + Breed.Group..1.', 'bray')## PERMANOVA Model: bray~Household + Breed.Group..1.; Strata: None## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 45 11.0303 0.69677 2.3174 9.999e-05 ***
## Breed.Group..1. 4 0.6750 0.04264 1.5954 0.0143 *
## Residual 39 4.1252 0.26059
## Total 88 15.8306 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1permdisp(ps.rare.prop, 'Breed.Group..1.', 'bray')## PERMDISP Model: bray~Breed.Group..1.
## missing observations due to 'group' removed##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 6 0.08392 0.013987 0.8302 10000 0.5411
## Residuals 82 1.38162 0.016849Weighted-Unifrac distance
plot_ord(ps.rare.prop, 'Breed.Group..1.','MDS','wunifrac')## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.rare.prop, 'Breed.Group..1.','NMDS','wunifrac')## Run 0 stress 0.1672646
## Run 1 stress 0.1682008
## Run 2 stress 0.1697386
## Run 3 stress 0.1667988
## ... New best solution
## ... Procrustes: rmse 0.06250708 max resid 0.1959744
## Run 4 stress 0.1700845
## Run 5 stress 0.1707173
## Run 6 stress 0.1711679
## Run 7 stress 0.174458
## Run 8 stress 0.1698938
## Run 9 stress 0.1651893
## ... New best solution
## ... Procrustes: rmse 0.06289663 max resid 0.2473875
## Run 10 stress 0.1707199
## Run 11 stress 0.1681964
## Run 12 stress 0.1766765
## Run 13 stress 0.1696632
## Run 14 stress 0.1684377
## Run 15 stress 0.1749852
## Run 16 stress 0.1769908
## Run 17 stress 0.1729453
## Run 18 stress 0.1654932
## ... Procrustes: rmse 0.066596 max resid 0.2438852
## Run 19 stress 0.1674721
## Run 20 stress 0.1698115
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 2: no. of iterations >= maxit
## 18: stress ratio > sratmax## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
permanova(ps.rare.prop, 'Household + Breed.Group..1.', 'wunifrac')## PERMANOVA Model: wunifrac~Household + Breed.Group..1.; Strata: None## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 45 1.72070 0.68078 2.2612 9.999e-05 ***
## Breed.Group..1. 4 0.14736 0.05830 2.1786 0.0041 **
## Residual 39 0.65950 0.26092
## Total 88 2.52756 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1permdisp(ps.rare.prop, 'Breed.Group..1.', 'wunifrac')## PERMDISP Model: wunifrac~Breed.Group..1.
## missing observations due to 'group' removed## Warning in betadisper(dist.matrix, group = df[, group]): some squared distances
## are negative and changed to zero##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 6 0.00979 0.0016310 0.341 10000 0.9156
## Residuals 82 0.39216 0.0047825Drug Effect
Alpha Diversity
sam.drug <- sam %>% filter(Epileptic.or.Control == 'Epileptic')
ggplot(sam.drug, aes(x = Pheno.Y.N, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
anova_test(Shannon~Pheno.Y.N, data = sam.drug, type = 1)## ANOVA Table (type I tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Pheno.Y.N 1 47 0.972 0.329 0.02Beta Diversity
Bray-Curtis distance
ps.drug <- ps.rare.prop %>% subset_samples(Epileptic.or.Control == 'Epileptic')
plot_ord(ps.drug, 'Pheno.Y.N','MDS','bray')
plot_ord(ps.drug, 'Pheno.Y.N','NMDS','bray')## Run 0 stress 0.2023653
## Run 1 stress 0.204988
## Run 2 stress 0.2063177
## Run 3 stress 0.200578
## ... New best solution
## ... Procrustes: rmse 0.02550675 max resid 0.1203173
## Run 4 stress 0.1998504
## ... New best solution
## ... Procrustes: rmse 0.06159441 max resid 0.3460622
## Run 5 stress 0.2067472
## Run 6 stress 0.2200867
## Run 7 stress 0.2018181
## Run 8 stress 0.1991884
## ... New best solution
## ... Procrustes: rmse 0.06579866 max resid 0.3503428
## Run 9 stress 0.2042939
## Run 10 stress 0.2048211
## Run 11 stress 0.2066389
## Run 12 stress 0.2227636
## Run 13 stress 0.2091282
## Run 14 stress 0.210517
## Run 15 stress 0.2055955
## Run 16 stress 0.2369509
## Run 17 stress 0.1998499
## Run 18 stress 0.210047
## Run 19 stress 0.2063561
## Run 20 stress 0.2064086
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 4: no. of iterations >= maxit
## 16: stress ratio > sratmax
permanova(ps.drug, 'Pheno.Y.N', 'bray')## PERMANOVA Model: bray~Pheno.Y.N; Strata: None## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Pheno.Y.N 1 0.2174 0.02545 1.2275 0.1978
## Residual 47 8.3243 0.97455
## Total 48 8.5417 1.00000permdisp(ps.drug, 'Pheno.Y.N', 'bray')## PERMDISP Model: bray~Pheno.Y.N##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.02021 0.020215 1.5701 10000 0.2261
## Residuals 47 0.60514 0.012875Weighted-Unifrac distance
plot_ord(ps.drug, 'Pheno.Y.N','MDS','wunifrac')
plot_ord(ps.drug, 'Pheno.Y.N','NMDS','wunifrac')## Run 0 stress 0.1777287
## Run 1 stress 0.1870529
## Run 2 stress 0.1816831
## Run 3 stress 0.1886616
## Run 4 stress 0.1881457
## Run 5 stress 0.1831296
## Run 6 stress 0.1902225
## Run 7 stress 0.177609
## ... New best solution
## ... Procrustes: rmse 0.04874602 max resid 0.1841465
## Run 8 stress 0.1950317
## Run 9 stress 0.180222
## Run 10 stress 0.1770123
## ... New best solution
## ... Procrustes: rmse 0.02868657 max resid 0.1280371
## Run 11 stress 0.1903529
## Run 12 stress 0.1798196
## Run 13 stress 0.1852662
## Run 14 stress 0.1959345
## Run 15 stress 0.1991224
## Run 16 stress 0.1795676
## Run 17 stress 0.1830592
## Run 18 stress 0.189317
## Run 19 stress 0.2013178
## Run 20 stress 0.1895554
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
permanova(ps.drug, 'Pheno.Y.N', 'wunifrac')## PERMANOVA Model: wunifrac~Pheno.Y.N; Strata: None## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Pheno.Y.N 1 0.02658 0.02163 1.0391 0.3848
## Residual 47 1.20210 0.97837
## Total 48 1.22868 1.00000permdisp(ps.drug, 'Pheno.Y.N', 'wunifrac')## PERMDISP Model: wunifrac~Pheno.Y.N##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.011062 0.0110619 3.6224 10000 0.06309 .
## Residuals 47 0.143525 0.0030537
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sex Effect
prop.test(xtabs(~Household+Sex, data = sam))## Warning in prop.test(xtabs(~Household + Sex, data = sam)): Chi-squared
## approximation may be incorrect##
## 49-sample test for equality of proportions without continuity
## correction
##
## data: xtabs(~Household + Sex, data = sam)
## X-squared = 52.464, df = 48, p-value = 0.3051
## alternative hypothesis: two.sided
## sample estimates:
## prop 1 prop 2 prop 3 prop 4 prop 5 prop 6 prop 7 prop 8 prop 9 prop 10
## 0.5 0.5 0.0 1.0 0.0 1.0 1.0 0.5 1.0 1.0
## prop 11 prop 12 prop 13 prop 14 prop 15 prop 16 prop 17 prop 18 prop 19 prop 20
## 0.5 1.0 1.0 1.0 0.5 0.5 1.0 1.0 0.0 0.5
## prop 21 prop 22 prop 23 prop 24 prop 25 prop 26 prop 27 prop 28 prop 29 prop 30
## 0.0 0.5 1.0 0.5 0.5 0.0 0.0 0.0 0.5 0.5
## prop 31 prop 32 prop 33 prop 34 prop 35 prop 36 prop 37 prop 38 prop 39 prop 40
## 1.0 0.5 1.0 1.0 0.0 0.5 0.5 0.5 0.5 0.5
## prop 41 prop 42 prop 43 prop 44 prop 45 prop 46 prop 47 prop 48 prop 49
## 0.0 1.0 1.0 0.5 0.5 1.0 0.5 1.0 0.5Alpha Diversity
ggplot(sam, aes(x = Sex, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
anova_test(Shannon~Household+Sex, data = sam, type = 1)## ANOVA Table (type I tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Household 48 48 1.627 0.048 * 0.619
## 2 Sex 1 48 2.087 0.155 0.042Beta Diversity
Bray-Curtis distance
plot_ord(ps.rare.prop, 'Sex','MDS','bray')
plot_ord(ps.rare.prop, 'Sex','NMDS','bray')## Run 0 stress 0.2088385
## Run 1 stress 0.2031576
## ... New best solution
## ... Procrustes: rmse 0.04283347 max resid 0.2838237
## Run 2 stress 0.2244714
## Run 3 stress 0.2041024
## Run 4 stress 0.2042066
## Run 5 stress 0.2039633
## Run 6 stress 0.2367018
## Run 7 stress 0.203949
## Run 8 stress 0.2162262
## Run 9 stress 0.203362
## ... Procrustes: rmse 0.01065205 max resid 0.07864169
## Run 10 stress 0.2038468
## Run 11 stress 0.2046384
## Run 12 stress 0.2045853
## Run 13 stress 0.2143459
## Run 14 stress 0.2043244
## Run 15 stress 0.2204904
## Run 16 stress 0.2031447
## ... New best solution
## ... Procrustes: rmse 0.01635909 max resid 0.08530236
## Run 17 stress 0.2095042
## Run 18 stress 0.2028724
## ... New best solution
## ... Procrustes: rmse 0.01143336 max resid 0.07831489
## Run 19 stress 0.2088665
## Run 20 stress 0.2232441
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 6: no. of iterations >= maxit
## 14: stress ratio > sratmax
permanova(ps.rare.prop, 'Household+Sex', 'bray')## PERMANOVA Model: bray~Household+Sex; Strata: None## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 12.2401 0.68965 2.2442 9.999e-05 ***
## Sex 1 0.0541 0.00305 0.4758 0.9874
## Residual 48 5.4541 0.30730
## Total 97 17.7482 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1permdisp(ps.rare.prop, 'Sex', 'bray')## PERMDISP Model: bray~Sex##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.01816 0.018164 1.2295 10000 0.2755
## Residuals 96 1.41828 0.014774Weighted-Unifrac distance
plot_ord(ps.rare.prop, 'Sex','MDS','wunifrac')
plot_ord(ps.rare.prop, 'Sex','NMDS','wunifrac')## Run 0 stress 0.1672646
## Run 1 stress 0.1719031
## Run 2 stress 0.1744687
## Run 3 stress 0.1679213
## Run 4 stress 0.1749359
## Run 5 stress 0.1687601
## Run 6 stress 0.1731133
## Run 7 stress 0.1759072
## Run 8 stress 0.1684304
## Run 9 stress 0.1724602
## Run 10 stress 0.1770518
## Run 11 stress 0.1654213
## ... New best solution
## ... Procrustes: rmse 0.06227521 max resid 0.2380526
## Run 12 stress 0.1676802
## Run 13 stress 0.1834458
## Run 14 stress 0.1753542
## Run 15 stress 0.1846202
## Run 16 stress 0.1700225
## Run 17 stress 0.1721061
## Run 18 stress 0.1764543
## Run 19 stress 0.167466
## Run 20 stress 0.1681448
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
permanova(ps.rare.prop, 'Household+Sex', 'wunifrac')## PERMANOVA Model: wunifrac~Household+Sex; Strata: None## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 1.90061 0.67646 2.0995 9.999e-05 ***
## Sex 1 0.00377 0.00134 0.1999 0.9977
## Residual 48 0.90527 0.32220
## Total 97 2.80965 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1permdisp(ps.rare.prop, 'Sex', 'wunifrac')## PERMDISP Model: wunifrac~Sex##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00156 0.0015650 0.391 10000 0.5351
## Residuals 96 0.38426 0.0040027Age Effect
Alpha Diversity
sam[which(is.na(sam$Age..months.)),'Household'] ## [1] "9" "24"# remove households that have dog with unspecific age
sam.Age <- sam %>% filter(!(Household %in% c('9', '24')))
ggplot(sam.Age,) +
geom_line(aes(x = as.numeric(Household), y = Age..months., group = Household)) +
geom_point(aes(x = as.numeric(Household), y = Age..months., group = Household, colour = Epileptic.or.Control)) +
xlab('Household') + ylab('Age in month')
anova_test(Shannon~Household+Age..months., data = sam.Age, type = 1)## ANOVA Table (type I tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Household 46 46 1.337 0.164 0.572000
## 2 Age..months. 1 46 0.024 0.877 0.000528Beta Diversity
Bray-Curtis distance
ps.age <- ps.rare.prop %>% subset_samples(!(Household %in% c('9', '24')))
permanova(ps.age, 'Age..months.', 'bray', strata = 'Household')## PERMANOVA Model: bray~Age..months.; Strata: Household## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Age..months. 1 0.274 0.01599 1.4951 0.296
## Residual 92 16.857 0.98401
## Total 93 17.131 1.00000Weighted-Unifrac distance
permanova(ps.age, 'Age..months.', 'wunifrac')## PERMANOVA Model: wunifrac~Age..months.; Strata: None## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, by = "term", parallel = core, na.action = na.omit, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Age..months. 1 0.03874 0.01428 1.3324 0.1956
## Residual 92 2.67474 0.98572
## Total 93 2.71348 1.00000sessioninfo::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.4.1 (2024-06-14)
## os macOS 15.1.1
## system aarch64, darwin20
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/New_York
## date 2024-12-29
## pandoc 3.5 @ /Users/yixuanyang/miniforge3/bin/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-8 2024-09-12 [1] CRAN (R 4.4.1)
## ade4 1.7-22 2023-02-06 [1] CRAN (R 4.4.0)
## ape 5.8-1 2024-12-16 [1] CRAN (R 4.4.1)
## backports 1.5.0 2024-05-23 [1] CRAN (R 4.4.0)
## Biobase 2.66.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## BiocGenerics 0.52.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## biomformat 1.34.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## Biostrings 2.74.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## broom 1.0.7 2024-09-26 [1] CRAN (R 4.4.1)
## bslib 0.8.0 2024-07-29 [1] CRAN (R 4.4.0)
## cachem 1.1.0 2024-05-16 [1] CRAN (R 4.4.0)
## car 3.1-3 2024-09-27 [1] CRAN (R 4.4.1)
## carData 3.0-5 2022-01-06 [1] CRAN (R 4.4.0)
## cli 3.6.3 2024-06-21 [1] CRAN (R 4.4.0)
## cluster 2.1.8 2024-12-11 [1] CRAN (R 4.4.1)
## codetools 0.2-20 2024-03-31 [1] CRAN (R 4.4.1)
## colorspace 2.1-1 2024-07-26 [1] CRAN (R 4.4.0)
## crayon 1.5.3 2024-06-20 [1] CRAN (R 4.4.0)
## data.table 1.16.4 2024-12-06 [1] CRAN (R 4.4.1)
## digest 0.6.37 2024-08-19 [1] CRAN (R 4.4.1)
## dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.4.0)
## evaluate 1.0.1 2024-10-10 [1] CRAN (R 4.4.1)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.4.0)
## fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.4.0)
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.4.0)
## foreach 1.5.2 2022-02-02 [1] CRAN (R 4.4.0)
## Formula 1.2-5 2023-02-24 [1] CRAN (R 4.4.0)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.4.0)
## GenomeInfoDb 1.42.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## GenomeInfoDbData 1.2.13 2024-10-05 [1] Bioconductor
## ggplot2 * 3.5.1 2024-04-23 [1] CRAN (R 4.4.0)
## glue 1.8.0 2024-09-30 [1] CRAN (R 4.4.1)
## gtable 0.3.6 2024-10-25 [1] CRAN (R 4.4.1)
## here * 1.0.1 2020-12-13 [1] CRAN (R 4.4.0)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.4.0)
## htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.4.0)
## httr 1.4.7 2023-08-15 [1] CRAN (R 4.4.0)
## igraph 2.1.2 2024-12-07 [1] CRAN (R 4.4.1)
## IRanges 2.40.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## iterators 1.0.14 2022-02-05 [1] CRAN (R 4.4.0)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.4.0)
## jsonlite 1.8.9 2024-09-20 [1] CRAN (R 4.4.1)
## knitr 1.49 2024-11-08 [1] CRAN (R 4.4.1)
## labeling 0.4.3 2023-08-29 [1] CRAN (R 4.4.0)
## lattice * 0.22-6 2024-03-20 [1] CRAN (R 4.4.1)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.4.0)
## lubridate * 1.9.4 2024-12-08 [1] CRAN (R 4.4.1)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.4.0)
## MASS 7.3-61 2024-06-13 [1] CRAN (R 4.4.0)
## Matrix 1.7-1 2024-10-18 [1] CRAN (R 4.4.1)
## mgcv 1.9-1 2023-12-21 [1] CRAN (R 4.4.1)
## multtest 2.62.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## munsell 0.5.1 2024-04-01 [1] CRAN (R 4.4.0)
## nlme 3.1-166 2024-08-14 [1] CRAN (R 4.4.0)
## permute * 0.9-7 2022-01-27 [1] CRAN (R 4.4.0)
## phyloseq * 1.50.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## pillar 1.10.0 2024-12-17 [1] CRAN (R 4.4.1)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.4.0)
## plyr 1.8.9 2023-10-02 [1] CRAN (R 4.4.0)
## purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.4.0)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.4.0)
## Rcpp 1.0.13-1 2024-11-02 [1] CRAN (R 4.4.1)
## readr * 2.1.5 2024-01-10 [1] CRAN (R 4.4.0)
## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.4.0)
## rhdf5 2.49.0 2024-05-04 [1] Bioconductor 3.20 (R 4.4.0)
## rhdf5filters 1.17.0 2024-05-04 [1] Bioconductor 3.20 (R 4.4.0)
## Rhdf5lib 1.27.0 2024-05-04 [1] Bioconductor 3.20 (R 4.4.0)
## rlang 1.1.4 2024-06-04 [1] CRAN (R 4.4.0)
## rmarkdown 2.29 2024-11-04 [1] CRAN (R 4.4.1)
## rprojroot 2.0.4 2023-11-05 [1] CRAN (R 4.4.0)
## rstatix * 0.7.2 2023-02-01 [1] CRAN (R 4.4.0)
## rstudioapi 0.17.1 2024-10-22 [1] CRAN (R 4.4.1)
## S4Vectors 0.44.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## sass 0.4.9 2024-03-15 [1] CRAN (R 4.4.0)
## scales 1.3.0 2023-11-28 [1] CRAN (R 4.4.0)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.4.0)
## stringi 1.8.4 2024-05-06 [1] CRAN (R 4.4.0)
## stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.4.0)
## survival 3.8-3 2024-12-17 [1] CRAN (R 4.4.1)
## tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.4.0)
## tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.4.0)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.4.0)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.4.0)
## timechange 0.3.0 2024-01-18 [1] CRAN (R 4.4.0)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.4.0)
## UCSC.utils 1.2.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.4.0)
## vegan * 2.6-8 2024-08-28 [1] CRAN (R 4.4.1)
## withr 3.0.2 2024-10-28 [1] CRAN (R 4.4.1)
## xfun 0.49 2024-10-31 [1] CRAN (R 4.4.1)
## XVector 0.46.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## yaml 2.3.10 2024-07-26 [1] CRAN (R 4.4.0)
## zlibbioc 1.52.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##
## [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
##
## ──────────────────────────────────────────────────────────────────────────────